Impact ionization rates for Si, GaAs, InAs, ZnS, and GaN in the GW approximation 
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We present first-principles calculations of the impact ionization rate (IIR) in the GW approxi- 
mation (GWA) for semiconductors. The IIR is calculated from the quasiparticle (QP) width in the 
GWA, since it can be identified as the decay rate of a QP into lower energy QP plus an independent 
electron-hole pair. The quasiparticle self-consistent GW method was used to generate the nonin- 
teracting hamiltonian the GWA requires as input. Small empirical corrections were added so as to 
reproduce experimental band gaps. Our results are in reasonable agreement with previous work, 
though we observe some discrepancy. In particular we find high IIR at low energy in the narrow 
gap semiconductor InAs. 
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I. INTRODUCTION 

The electron-initiated impact ionization is a fundamen- 
tal process in semiconductors where a high energy elec- 
tron decays into an another low-energy electron together 
with an electron- hole pair [1]. The impact ionization rate 
(IIR), which originates from the coulomb interaction be- 
tween electrons, is a critical factor affecting transport 
under high electric field, as described by the Boltzmann 
transport equation (BTE). It is important in narrow gap 
semiconductors, especially for ultrasmall devices. Impact 
ionization is also used in avalanche photodiodes, and to 
supply electron- hole pairs for electroluminescence. Re- 
cently it has stimulated interest as a mechanism to im- 
prove efficiency in photovoltaic devices [2[. 

The IIR has been calculated with empirical pseudopo- 
tentials (EPP) in order to include realistic energy bands 
B PL Si B 3- Sano and Yoshii calculated the IIR for 
Si [4], [5[ and obtained reasonable agreement with ex- 
perimental data. They also studied other materials @, 
treating the transition matrix element M as a parameter 
(constant matrix approximation). Jung et al. [7] used 
an EPP to calculate the IIR in GaAs. They calculated 
M including explicit calculation of the dielectric function 
e(q, cj), rather than assuming a model form. 

Recently, two groups have calculated the IRR using the 
density-functional formalism to generate the one-body 
eigenfunctions and energy bands. Because the standard 
local density approximation (LDA) underestimates semi- 
conductor bandgaps while the IIR is very sensitive to 
this quantity, the standard LDA is not suitable. Pi- 
cozzi et al. used a screened-exchange generalization [ij 
of the LDA l9UlOL and Kuligk et al. employed the exact 
exchange jl ill 12| formalism [13[. Both groups used model 
dielectric functions for the dynamically screened coulomb 
interaction W(r, r', uo). 

Here we will present (nearly) ab initio calculations of 
the IIR without model assumptions. First, our noninter- 
acting hamiltonian H° is generated within quasiparticle 
self-consistent GW (QS GW) formalism. We have shown 
that QS G W works very well for wide range of materials 
[H El EE E3, O- Because the IIR is highly sensi- 



tive to the bandgap, we add a small empirical scaling 
of the exchange-correlation potential so as to reproduce 
the experimental fundamental gap Eg- Corrections for 
semiconductors are small and systematic as shown below. 
Second, W is calculated from the QS GW noninteracting 
hamiltonian. The IIR is identified with the decay rate 
(or linewidth) of the quasiparticle (QP), which is calcu- 
lated from the imaginary part of the self-energy, as we 
describe below. Our method thus contains only one pa- 
rameter, to correct the band gap. As we have shown [19j, 
this parameter is small and is approximately independent 
of material. In principle our method can predict the IIR 
in unknown systems, and also for inhomogeneous systems 
such as grain boundaries, quantum dots, or impurities, 
where the IIR should be strongly enhanced because mo- 
mentum conservation is much more easily satisfied. Thus 
the present ab initio method should be superior to prior 
approaches. Applications to such systems will be useful 
in devices that need to suppress or enhance electron-hole 
pair-generation from impact ionization. 

After a theoretical discussion, we present some results. 
They are in reasonable agreement with previous calcula- 
tions, except for InAs where IIR is calculated to be much 
higher than what Sano and Yoshii found @. 



II. METHOD 

The first step is to determine a good one-body Hamil- 
tonian H° which describes QPs. We obtain H° from 
QS GW calculations 0, [H, Il6j . As we explain in Sec lIIIi 
we follow Ref.[19], and modify H° by a simple empirical 
scaling (a-correction) to ensure the fundamental gap re- 
produces experiment. From this modified H° we obtain a 
set of eigenvalues {^kn} and eigenfunctions {^kn} 7 which 
are used to calculate the self-energy E(r, r', uo) within the 
GWA, Z = iGxW. The inverse of the QP lifetime r^ 1 
is obtained from the imaginary part of E as 
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where ImE kn = (^kn|I m ^(^kn)|^kn) = e.g, see Eq.(l) in Ref.Q. It is written as 
/d 3 r/d 3 r / ^* n (r)ImS(r,r / ,e kn )^ kn (r / ). By ImE(s kn ) 
we mean the anti-hermitian part. Z kn is the wave 
function renormalization factor to represent the QP 
weight, k denotes the wave vector in the first Brillouin 
zone (BZ), and n the band index. The expression EqU] 
for is derived in Appendix B of Ref.[20]. ImE is 
obtained from the the imaginary part of the convolution 
of G and W. For an unoccupied state kn, it is 
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where states k'n' are restricted to those for which ef < 
£ k / n / < e kn . is calculated in the random phase ap- 
proximation (RPA) as 



W = v + = (1 — ^Xo) 1 v, 



(3) 



where v is the coulomb interaction; x is the full polariza- 
tion function in the RPA, and xo is the non-interacting 
polarization function. With Eqs. (|l|2l3p , r^ 1 is calculated 
from H° in principle. In the Lehmann representation, x 
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where |ra) denotes the eigenstates (intermediate states) 
with excitation energy uj m relative to the ground state 
|0). Here n(r) is the density operator. 

In the RPA, \m) are the eigenfunctions of a two-body 
(one electron and one hole) eigenvalue problem in the 
RPA. In simple cases such as the homogeneous electron 
gas, \m) for high uo m are identified as plasmons; \m) for 
low com are as independent motions of an electron and 
a hole. Thus r^ 1 for low energy electrons calculated in 
GWA can be identified as the transition probability to 
such states for the independent motion of an electron and 
a hole together with an electron; that is, we identify r^ 1 
as the IIR. 

There are some questionable points for the identifica- 
tion. It might be not so easy in some cases to iden- 
tify a state \m) as such a independent motion because 
the electron-hole pair can be hybridized with plasmons. 
However, such hybridization is sufficiently small for the 
simple semiconductors treated here, because plasmons 
appear only at high energies as u m > 1 Ry. Another 
problem is that the final state consisting of two electrons 
and one hole is not symmetrized for the electrons in the 
GWA. Thus Fermi statistics are not satisfied. Below we 
discuss how much error it causes. 

Our formula Eq. (pQ) for r^ 1 is different from the cus- 
tomary expression found in the literature 0, 0, H, d, 5 



where |M| 2 

eludes both direct and exchange processes. The sum over 
k',ki,k2 is restricted to satisfy k = k 7 + ki — k2. The 
matrix element Md for the direct process is 
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Me for the exchange process is the same as Md, ex- 
cept that the two electrons in final states (k'n' kini) 
are exchanged. Eq. ([5]) can be derived in time-dependent 
perturbation theory, where the final states consists of two 
electrons and one hole. This is based on the physical pic- 
ture that W causes transitions between the Fock states 
made of QPs. However, the final states made of the three 
QPs are interacting each other. Thus such a picture do 
not necessarily well-defined. This is related to a funda- 
mental problem about how to mimic the quantum theory 
by the BTE. Definition of the IIR suitable for the BTE is 
somehow ambiguous. The difference between Eq. (pQ) and 
Eq. ([5]) is related to the ambiguity. One is not necessarily 
better than the other. 

To compare Eq. (pQ) with Eq. ([5]), let us assume that 
Imxo is small enough. Then we have 



ImW ~ Wr Imxo Wr 



(7) 



from Eq. (|3]), where Wr = (1 — v~Rexo)~ 1 v. Rexo denotes 
the hermitian (real) part of xo- If we apply Eq. (JJ) to 
Eq. ([2j), Eq. (pQ) is reduced to an expression similar to 
Eq. ©: 
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where Md is defined in Eq. (j6]) but with Wr instead of 
W. Through Eq. (jHj) we can elucidate the differences 
between Eq. (pQ) and Eq. ([5|) as follows: 

(a) Eq. J5J) (and thus Eq. (pQ)) contain the Z factor. 
This is because Eq. ([5]) was derived without tak- 
ing into account the modification of QPs by the 
coulomb interaction. Typically Z kn is ~ 0.8. 

(b) Eq. ((5]) and Eq. (pQ) do not include Me contribu- 
tions. In the extreme case when Me = ^Md , 
|M| 2 = 0.75 x |M D | 2 . This occurs in the Hubbard 
model when W is a point interaction; the Feyn- 
man diagrams for Md and Me become the same 
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TABLE I: Eigenvalues of semiconductors relative to the va- 
lence band maximum at F. The QSGW column depict re- 
sults without spin-orbit coupling; values in parentheses in- 
clude spin-orbit coupling. Column QSGWa shows values af- 
ter scaling defined in Eq. J9}. a is chosen so that the QSGWa 
potential (without spin-orbit coupling) reproduces the exper- 
imental minimum band gap at room temperature. 
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Experimental data at room temperature, taken from Ref. |2l). 
c Empirical pseudopotential data are taken from Ref. |2l). The 
experimental direct gap is ~0.4eV. 



except for their sign. Theoretically, including Me is 
advantageous because it symmetrizes the two elec- 
trons in the final state (though only for Imxo in 
the linear response regime). Fermi statistics are not 
perfectly satisfied because not all the exchange-pair 
diagrams are included. Omitting the exchange con- 
tribution reduces the IIR by a factor 0.75 at most, 
as explained above. 



(c) Eq. (|8j) contains only the real part of W, in contrast 
to Eq. ([5]) . The difference originates from higher or- 
der contributions to Imxo- Moreover, when Eq. ([7j) 
is not satisfied there are further higher-order con- 
tributions to Imxo- 



We may have to pay attention to these differences. For 
small Imxo? ( a ) an d (b) predominate, and the difference 
between Eq. (pQ) and the Kane formula Eq. ([5]) should be 
a factor in the range 0.5 to 1. However, this difference is 
relatively minor on the log scale in the Figure. 



III. RESULTS 

Here we treat Si, GaAs, InAs, zincblende ZnS, and 
wurtzite GaN. For each material we calculate a self- 
consistent noninteracting hamiltonian H° through the 
QSGW formalism. Spin-orbit coupling is neglected, fol- 
lowing prior work [6|, lZ[ . Table [J shows calculated values 
at high-symmetry points, compared with available exper- 
imental data. As we and others have noted [H, [l6|, [Hf, 
the QSGW gap is systematically overestimated because 
the RPA underestimates the screening. (Also the GaAs 
calculation used a smaller basis what was reported in 
[IH, resulting in an additional overestimate of ~0.05 eV.) 
To compare the QSGW results to experiment, we must 
take into account other contributions: spin-orbit cou- 
pling, zero-point motionj23[, and finite temperature all 
reduce the gap slightly [24[. 

To obtain the most reliable IIR, we slightly modify H° 
to reproduce the experimental gap at room temperature 
without including these contributions explicitly. To do 
this, we add an empirical scaling ("QSGWa" correction) 
following the procedure used in Ref. [191]. We scale the 
one-body Hamiltonian as follows: 

H a = H° + (1 - a)(£ - K L C DA ) (9) 

where E is the static version of the self-energy (see 
Eq.(10) in Ref. [16]). Table [J shows numerical values both 
with and without the scaling. As we showed in Ref.pjjj], 
effective masses are also well reproduced. Thus we can 
set up a satisfactory H a with a single parameter a. This 
procedure is reasonable because the uncorrected gaps are 
already close to experiment and 1 — a is not large. The 
materials-dependence of a shown in Table [J originates 
largely from the dependence of SO coupling and finite 
temperature on material. If these were taken into account 
by improving H° explicitly, a universal choice of a ~ 0.8 
would reproduce the experimental gaps in the Table to 
within ~0.1eV. (Alternatively, adopting the present pro- 
cedure with a universal a ~ 0.75 accomplishes much the 
same thing.) Table [J shows that the experimental en- 
ergy dispersions are also well reproduced where they are 
well known (Si and GaAs). This systematic tendency is 
found for many other materials, including ZnO, CU2O, 
NiO and MnO [MEl], and GdN0. It implies that the 
QSGWa procedure is broadly applicable with compara- 
ble accuracy to many environments, e.g, to InAs/GaAs 
grain boundaries. The QSGWa energy bands are shown 
in FigHJ 

Given H a , we perform a one-shot GWA calculation 
using the method detailed in Ref.[l6j], and calculate r^ 1 
from Eq. (pQ) . To reduce the computational time we trun- 
cate the product basis for each atomic site to I < 1. This 
limits the degrees of freedom for the local-field correc- 
tion in the dielectric function. However, we checked that 
this little affects the results. To obtain Z\^ n in Eq. (pQ), we 
need to calculate the derivative of the self-energy dRe ^}^ 
at £kn, though Z contributes a relatively unimportant 
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FIG. 1: Energy bands calculated by the QSGWa method (Eq. ©), with a chosen to reproduce the experimental band gap 
Eg- In wurtzite GaN, i?G=3.44eV [21], and a = 0.79. Data for other compounds can be found in Table HI 



factor ~ 0.8. The main computational cost of the IIR 
calculation comes from the sum of the pole weights on 
the real axis; see Eq. (58) in Ref.[16[. This corresponds 
to the convolution of ImG and ImW, Eq. after ImW 
is obtained from integration by the tetrahedron method 
pi| . Fig 12] shows our results for r^ 1 . The x-axis denotes 
the initial electron energy e- in measured from the bottom 
of the conduction band; e- in > Eg is a hard threshold 
below which IIR is zero. The present results, depicted 
by large plus signs, are superposed on results taken from 
previous work. 

The IIR has a typical feature as already shown in Fig.l 
m Ref.@, that is, the IIR as function of Si n are widely 
scattered at low S{ n because of the limited number of 
transitions that conserve energy and momentum. The 
scatter diminishes at high energy because of an averag- 
ing effect which smears the anisotropy in the Brillouin 
zone as discussed in [5[. Our results for Si and GaAs 
correspond rather well to previous work. Details for the 
IIR are already well analyzed 0, i, i, 0, S EO, . 

Turning to ZnS and GaN, we superpose our results 
on those presented by Kuligk et al. in Figs. 9 and 10 of 



Ref. pjl , which include exact exchange (EXX) (solid sym- 
bols) and EPP results from Ref. [10] (open symbols). The 
EPP and the present calculations appear mostly similar 
apart from an approximately constant factor; however 
the EXX results show rather different behavior, particu- 
larly in GaN. This is likely because the EPP and QS GWa 
energy bands are quite similar to each other, but they are 
quite different from the EXX case (see Figs 2 and 3 in 
Ref.[H). 

A large discrepancy with EPP is seen only in InAs. 
Our data is superposed on the calculations by Sano and 
Yoshii [6] . We obtain high IIR at low initial electron ener- 
gies S[ n ^1 eV. Such high IIR comes from initial electrons 
near the conduction band minimum at the T point. Since 
the band gap and effective mass are small in InAs, there 
are states not far from T with energy Si n > Eg, which 
can generate an electron-hole pair. This occurs only for 
InAs in the cases studied, but generally occurs for narrow 
gap semiconductors. For the discrepancy with results of 
Sano and Yoshii may be due to their constant matrix el- 
ements approximation, which is not suitable for such a 
narrow gap material (see Ref.@ near Eq.(2)). 
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FIG. 2: (color online) Impact ionization rates t^ 1 as a function of the initial electron energy Si n (measured from the bottom of 
conduction band). The present GWA calculation is shown by large (red) plus signs. It is superposed on previous calculations: 
Si and InAs from Sano and Yoshii [6], GaAs from Jung et al. Q], ZnS and GaN from Kuligk et a l. Il3| ]. Open boxes in the 
ZnS and GaN data are EPP results from Picozzi et al. [10| : solid circles are exact exchange results [l3|. We used 500 k points 
in the 1st Brillouin zone for GaN, and 1728 points for the cubic compounds (regular mesh including the T point [16]). Owing 
to the limited number of k points, there are some numerical errors, e.g. a factor of order 2 when IIR < le+10). The error is 
not large enough to affect our conclusions. 
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In conclusion, we have calculated the IIR for several 
materials in the GWA, after a theoretical discussion of 
its application to the IIR. In principle, the method pre- 
sented here will be applicable even to inhomogeneous sys- 
tems such grain boundaries and quantum dots where we 
expect very strong IIR. The present calculations corre- 
spond reasonably well to prior work, with the exception 
of the narrow gap material In As. High IIR would be ex- 



pected universally in similar narrow gap materials such 
as GaSb, InSb, and InN. This indicates that careful con- 
sideration for the IIR might be required when we use 
such materials for devices. 
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